
typedef struct _Interval {
	float low;
	float up;
} Interval;

float sqr(float x) {
	return x * x;
}

// Compute N(c), the number of eigenvalues less than c
int getNumEigensLess(float c, int n, __global float *d, __global float *e, float pivmin) {
	
	float diff = d[0] - c;
	diff = select(diff, -pivmin, fabs(diff) < pivmin);
	int count = select(0, 1, diff <= 0);

	for (int j = 1; j < n; ++j) {
		diff = d[j] - c - sqr(e[j - 1]) / diff;
		diff = select(diff, -pivmin, fabs(diff) < pivmin);
		count += select(0, 1, diff <= 0);
	}
	return count;
}

// returns true if interval is narrow enough to stop bisection its
int isNarrow(Interval interval, float abstol, float reltol) {
	float width = fabs(interval.up - interval.low);
	float norm = max(fabs(interval.up), fabs(interval.low));
	float minWidth = max(abstol, reltol * norm);
	return (width < minWidth);
}

__kernel void recalculateEigenInterval(int n,
		__global float *d, __global float *e,
		__global Interval* inIntervals, __global Interval* outIntervals,
		__global int *inNumLess, __global int *outNumLess,		
		float abstol, float reltol, float pivmin)
{
	uint threadId = get_global_id(0);

	int input = 0;
	int eigenNumber = threadId + inNumLess[0];
	while (inNumLess[input + 1] <= eigenNumber) {
		++input;
	}

	Interval in = inIntervals[input];
	
	int numLessLow = inNumLess[input];
	int numLessUp = inNumLess[input + 1];
		
	if (isNarrow(inIntervals[input], max(abstol, pivmin), reltol)) {
		// interval is converged - exit
		outIntervals[threadId] = in;
		outNumLess[threadId + 1] = numLessUp;
		return;
	}
	
	int numEigens = numLessUp - numLessLow;
	int useBisection = (numEigens == 1);
	
	float width = fabs(in.up - in.low);
	float middle = (in.up + in.low) * 0.5f;	
	
	int divisionIndex = threadId + inNumLess[0] - numLessLow;	
	float divisionWidth = width / numEigens;

	float divisionLow = in.low + divisionWidth * divisionIndex;
	float divisionUp = divisionLow + divisionWidth;
	
	float newPoint = select(divisionUp, middle, useBisection);	
	int numLessNew = getNumEigensLess(newPoint, n, d, e, pivmin);
	numLessNew = min(max(numLessNew, numLessLow), numLessUp); // insure that N(w) is monotone 
	
	int useUpperHalf = (numLessNew <= numLessLow);
	float bisectionLow = select(in.low, middle, useUpperHalf);
	float bisectionUp = select(middle, in.up, useUpperHalf);
	
	outIntervals[threadId].low = select(divisionLow, bisectionLow, useBisection);
	outIntervals[threadId].up = select(divisionUp, bisectionUp, useBisection);		
	outNumLess[threadId + 1] = select(numLessNew, numLessUp, useBisection && useUpperHalf);
}


// XXX

__kernel void calculateNumEigensInInterval(__global int *numEigens,
	__global Interval* intervals,
	int n, __global float *d, __global float *e,
	float pivmin) 
{
	uint i = get_global_id(0);
	int numLessLower = getNumEigensLess(intervals[i].low, n, d, e, pivmin);
	int numLessUpper = getNumEigensLess(intervals[i].up, n, d, e, pivmin);
	numEigens[i] = max(numLessUpper - numLessLower, 0);
}

__kernel void recalculateInterval(
		int n, __global float *d, __global float *e,
		__global Interval *inIntervals, __global Interval* outIntervals,
		__global int* numEigens,
		float abstol, float reltol, float pivmin
)
{
	int threadId = get_global_id(0);

	int splitIndex = threadId;
	int inIndex = 0;
	while (splitIndex >= numEigens[inIndex]) {
		splitIndex -= numEigens[inIndex];
		++inIndex;
	}
	
	if (numEigens[inIndex] == 1) {
		float middle = 0.5f * (inIntervals[inIndex].up + inIntervals[inIndex].low);
		int num = getNumEigensLess(middle, n, d, e, pivmin);
		num -= getNumEigensLess(inIntervals[inIndex].low, n, d, e, pivmin);

		if (isNarrow(inIntervals[inIndex], abstol, reltol)) {
			outIntervals[threadId] = inIntervals[inIndex];
		} else if (num == 0) {
			// eigenvalue in upper half
			outIntervals[threadId].low = middle;
			outIntervals[threadId].up = inIntervals[inIndex].up;
		} else {
			// eigenvalue in lower half
			outIntervals[threadId].low = inIntervals[inIndex].low;
			outIntervals[threadId].up = middle;
		}
	} else {
		// more than 1 value in interval
		float width = (inIntervals[inIndex].up - inIntervals[inIndex].low) / numEigens[inIndex];
        outIntervals[threadId].low = inIntervals[inIndex].low + width * splitIndex;
        outIntervals[threadId].up = outIntervals[threadId].low + width;
	}
}

